Points

snxy <- data.frame(name = "Space Needle", x = -122.3493, y = 47.6205) 
# x:longitude, y:latitude
space_needle <- st_as_sf(snxy, coords = c("x", "y"), crs = 4326)

print(space_needle)
## Simple feature collection with 1 feature and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -122.3493 ymin: 47.6205 xmax: -122.3493 ymax: 47.6205
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##           name                  geometry
## 1 Space Needle POINT (-122.3493 47.6205)
st_coordinates(space_needle)
##           X       Y
## 1 -122.3493 47.6205
shxy <- data.frame(name = "Savery Hall", x = -122.3083, y = 47.6572)
savery_hall <- st_as_sf(shxy, coords = c("x", "y"), crs = 4326)

# rbind() to put two points in one data frame
pts <- rbind(space_needle, savery_hall)
print(pts)
## Simple feature collection with 2 features and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -122.3493 ymin: 47.6205 xmax: -122.3083 ymax: 47.6572
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##           name                  geometry
## 1 Space Needle POINT (-122.3493 47.6205)
## 2  Savery Hall POINT (-122.3083 47.6572)
plot(pts$geometry, axes = TRUE)

Lines

# create a linestring sf data frame 
lnstr <- st_sfc(st_linestring(st_coordinates(pts)), crs = 4326)
lnstr <- as_tibble(lnstr) %>% mutate(od = "Space Needle, Savery Hall")

plot(pts$geometry, axes = TRUE)
text(x = st_coordinates(pts), labels = pts$name)
plot(lnstr$geometry, col = 2, add = TRUE)

## Polygons

zooxy <- data.frame(name = "Woodland Park Zoo", x = -122.3543, y = 47.6685)
wp_zoo <- st_as_sf(zooxy, coords = c("x", "y"), crs = 4326)

# rbind() to put two points in one data frame
pts <- rbind(pts, wp_zoo)

(plygn <- st_sfc(st_polygon(list(st_coordinates(rbind(pts, space_needle)))), crs = 4326))
## Geometry set for 1 feature 
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -122.3543 ymin: 47.6205 xmax: -122.3083 ymax: 47.6685
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## POLYGON ((-122.3493 47.6205, -122.3083 47.6572,...
plot(plygn, col = "cyan", axes = TRUE)
plot(lnstr$geometry, col = 2, add = TRUE, lwd = 3)
plot(pts$geometry, add = TRUE, cex = 2)
text(x = st_coordinates(pts), labels = pts$name)

## Importing spatial data sets

# path to the data
mydatadir <- file.path("C:", "users", Sys.getenv("USERNAME"), "Documents","study",  "UW", "GIS_R","data")
zippolyfname <- file.path(mydatadir, "zip_poly.gdb")
# avoid reading over and over
if(!exists("zipcodes")){
    zipcodes <- st_read(dsn = zippolyfname, layer = "zip_poly", as_tibble = TRUE, geometry_column = "Shape")
}
## Reading layer `zip_poly' from data source `C:\Users\Yohan_Min\Documents\study\UW\GIS_R\data\zip_poly.gdb' using driver `OpenFileGDB'
## Simple feature collection with 30924 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -179.1473 ymin: 13.23419 xmax: 179.7785 ymax: 71.39048
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
# change the data frame's column names to lowercase
colnames(zipcodes) <- tolower(colnames(zipcodes))
# after renaming columns it is necessary to re-establish which column contains the geometry
st_geometry(zipcodes) <- "shape"
zip_wa <- zipcodes %>% filter(state == "WA")
head(zip_wa)
## Simple feature collection with 6 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -124.0031 ymin: 46.0401 xmax: -119.2583 ymax: 49.0014
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 9
##   zip_code po_name state population pop_sqmi   sqmi shape_length shape_area
##   <fct>    <fct>   <fct>      <int>    <dbl>  <dbl>        <dbl>      <dbl>
## 1 00072    Nation~ WA           -99     -99  1.22e3        4.20     0.386  
## 2 00073    Usdoe ~ WA           -99     -99  3.24e2        1.63     0.0984 
## 3 00074    Yakima~ WA           -99     -99  1.35e3        4.51     0.408  
## 4 00076    Okanog~ WA           -99     -99  7.31e2        3.83     0.232  
## 5 00195    Long I~ WA           -99     -99  8.37e0        0.536    0.00254
## 6 98001    Auburn  WA         35721    2003. 1.78e1        0.628    0.00549
## # ... with 1 more variable: shape <MULTIPOLYGON [°]>
plot(x = zip_wa$shape, axes = TRUE)

hospitals <- st_read(file.path(mydatadir, "medical_facilities/medical_facilities.shp"))
## Reading layer `medical_facilities' from data source `C:\Users\Yohan_Min\Documents\study\UW\GIS_R\data\medical_facilities\medical_facilities.shp' using driver `ESRI Shapefile'
## Simple feature collection with 154 features and 14 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 1254366 ymin: 78085.2 xmax: 1395044 ymax: 287325.1
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=47.5 +lat_2=48.73333333333333 +lat_0=47 +lon_0=-120.8333333333333 +x_0=500000.0000000001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
h2o <- st_read(file.path(mydatadir, "wtrbdy/wtrbdy.shp"))
## Reading layer `wtrbdy' from data source `C:\Users\Yohan_Min\Documents\study\UW\GIS_R\data\wtrbdy\wtrbdy.shp' using driver `ESRI Shapefile'
## Simple feature collection with 15838 features and 16 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 607316.1 ymin: -256676.3 xmax: 1617446 ymax: 765087.4
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=47.5 +lat_2=48.73333333333333 +lat_0=47 +lon_0=-120.8333333333333 +x_0=500000.0000000001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
# type = "n" not to plot, but sets xlim and ylim
plot(hospitals$geometry, type = "n", axes = TRUE)
# water
plot(h2o$geometry, col = "cyan", border = 0, add = TRUE)
# ZIP code areas
plot(zip_wa$shape, add = TRUE, col = 0, border = 1)
# hospital points
plot(hospitals$geometry, add = TRUE)
box()

# different coordinates 
# st_crs(hospitals)
# st_crs(zip_wa)

Exporting data

# st_write(obj = zip_wa, dsn = file.path(mydatadir, "zip_wa.shp"))
st_write(obj = pts, dsn = file.path(mydatadir, "r_gis.gpkg"), layer = "pts")
## Updating layer `pts' to data source `C:/users/Yohan_Min/Documents/study/UW/GIS_R/data/r_gis.gpkg' using driver `GPKG'
## Updating existing layer pts
## Writing 3 features with 1 fields and geometry type Point.
st_write(obj = lnstr, dsn = file.path(mydatadir, "r_gis.gpkg"), layer = "lnstr")
## Updating layer `lnstr' to data source `C:/users/Yohan_Min/Documents/study/UW/GIS_R/data/r_gis.gpkg' using driver `GPKG'
## Updating existing layer lnstr
## Writing 1 features with 1 fields and geometry type Line String.
st_write(obj = plygn, dsn = file.path(mydatadir, "r_gis.gpkg"), layer = "plygn")
## Updating layer `plygn' to data source `C:/users/Yohan_Min/Documents/study/UW/GIS_R/data/r_gis.gpkg' using driver `GPKG'
## Updating existing layer plygn
## Writing 1 features with 0 fields and geometry type Polygon.
st_write(obj = zip_wa, dsn = file.path(mydatadir, "r_gis.gpkg"), layer = "zip_wa")
## Updating layer `zip_wa' to data source `C:/users/Yohan_Min/Documents/study/UW/GIS_R/data/r_gis.gpkg' using driver `GPKG'
## Updating existing layer zip_wa
## Writing 547 features with 8 fields and geometry type Multi Polygon.

Coordinates

snxy <- data.frame(name = "Space Needle", x = -122.3493, y = 47.6205)
space_needle <- st_as_sf(snxy, coords = c("x", "y"))
st_crs(space_needle)
## Coordinate Reference System: NA
st_crs(space_needle) <- 4326
st_crs(space_needle)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.4-4, (SVN revision 833)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/Yohan_Min/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Yohan_Min/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.3-1
epsg <- make_EPSG()
utm10 <- epsg[grep("UTM.*10", epsg$note),]
kable(utm10) %>% kable_styling(bootstrap_options = "striped", full_width = F, position = "left")
code note prj4
1647 3157 # NAD83(CSRS) / UTM zone 10N +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
2207 3717 # NAD83(NSRS2007) / UTM zone 10N +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
2230 3740 # NAD83(HARN) / UTM zone 10N +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
3031 6339 # NAD83(2011) / UTM zone 10N +proj=utm +zone=10 +ellps=GRS80 +units=m +no_defs
4126 26710 # NAD27 / UTM zone 10N +proj=utm +zone=10 +datum=NAD27 +units=m +no_defs
4274 26910 # NAD83 / UTM zone 10N +proj=utm +zone=10 +datum=NAD83 +units=m +no_defs
4850 32210 # WGS 72 / UTM zone 10N +proj=utm +zone=10 +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.2263 +units=m +no_defs
4910 32310 # WGS 72 / UTM zone 10S +proj=utm +zone=10 +south +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.2263 +units=m +no_defs
4970 32410 # WGS 72BE / UTM zone 10N +proj=utm +zone=10 +ellps=WGS72 +towgs84=0,0,1.9,0,0,0.814,-0.38 +units=m +no_defs
5030 32510 # WGS 72BE / UTM zone 10S +proj=utm +zone=10 +south +ellps=WGS72 +towgs84=0,0,1.9,0,0,0.814,-0.38 +units=m +no_defs
5091 32610 # WGS 84 / UTM zone 10N +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
5159 32710 # WGS 84 / UTM zone 10S +proj=utm +zone=10 +south +datum=WGS84 +units=m +no_defs
5470 6653 # NAD83(CSRS) / UTM zone 10N + CGVD2013 height +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +vunits=m +no_defs 
kable(epsg %>%filter(code == 2927)) %>% kable_styling(bootstrap_options = "striped", full_width = F, position = "left")
code note prj4
2927 # NAD83(HARN) / Washington South (ftUS) +proj=lcc +lat_1=47.33333333333334 +lat_2=45.83333333333334 +lat_0=45.33333333333334 +lon_0=-120.5 +x_0=500000.0001016001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs 
(space_needle_utm10 <- space_needle %>% st_transform(26910))
## Simple feature collection with 1 feature and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 548894.1 ymin: 5274327 xmax: 548894.1 ymax: 5274327
## epsg (SRID):    26910
## proj4string:    +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
##           name                 geometry
## 1 Space Needle POINT (548894.1 5274327)
(space_needle <- space_needle %>% st_transform(26910))
## Simple feature collection with 1 feature and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 548894.1 ymin: 5274327 xmax: 548894.1 ymax: 5274327
## epsg (SRID):    26910
## proj4string:    +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
##           name                 geometry
## 1 Space Needle POINT (548894.1 5274327)
# st_crs(hospitals)
# st_crs(zip_wa)
# st_crs(h2o)

h2o <- h2o %>% st_transform(4326)
hospitals <- hospitals %>% st_transform(4326)

Geoprocessing

Buffering

# create the points
snxy <- data.frame(name = "Space Needle", x = -122.3493, y = 47.6205)
space_needle <- st_as_sf(snxy, coords = c("x", "y"), crs = 4326)
shxy <- data.frame(name = "Savery Hall", x = -122.3083, y = 47.6572)
savery_hall <- st_as_sf(shxy, coords = c("x", "y"), crs = 4326)
zooxy <- data.frame(name = "Woodland Park Zoo", x = -122.3543, y = 47.6685)
wp_zoo <- st_as_sf(zooxy, coords = c("x", "y"), crs = 4326)
pts <- rbind(space_needle, savery_hall, wp_zoo)

# make the buffer with inline transforms
pts_buf_1km <- pts %>% st_transform(26910) %>% st_buffer(dist = 1000) %>% st_transform(2926)
# crs to be converted to Cartesian projected coordinate for the same measure of distance. 

# write to the GPKG
mygpkg <- file.path(mydatadir, "r_gis.gpkg")
st_write(obj = pts_buf_1km, dsn = mygpkg, layer = "pts_buf_1km", quiet = TRUE, update = TRUE)
if(! exists("kctrans")){
    kctrans <- st_read(
        file.path(mydatadir, "Metro_Transportation_Network_TNET_in_King_County__trans_network_line.shp"),
        quiet = TRUE)
}

# freeways are KC_FCC_ID = F
kcfwy <- kctrans %>% filter(KC_FCC_ID == "F")

# buffer
kcfwy_buf_500ft <- kcfwy %>% st_transform(2926) %>% st_buffer(500)

# write to the GPKG
mygpkg <- file.path(mydatadir, "r_gis.gpkg")
st_write(obj = kcfwy_buf_500ft, dsn = mygpkg, layer = "kcfwy_buf_500ft", quiet = TRUE, update = TRUE)

Point-in-polygon

# if no API key,
# acs5_2018_bg <- st_read(dsn = file.path(mydatadir, "census.gpkg"), layer = "acs5_2018_bg", quiet = TRUE)
# st_crs(acs5_2018_bg) <- 2926

# cache data
options(tigris_use_cache = TRUE)
# where to store data
tigris_cache_dir <- mydatadir

# if you have your API key, enter it here rather than using the system environment variable
# myapikey <- "foobar"
myapikey <- Sys.getenv("CENSUS_API_KEY")
census_api_key(myapikey)
## To install your API key for use in future sessions, run this function with `install = TRUE`.
# get the data and project it to match the bus stops, also calculate the area
acs5_2018_bg <- get_acs(
    geography = "block group",
    variables = c(medfamincome="B19113_001"),
    state = "WA",
    county = "King",
    geometry = TRUE,
    moe = 95,
    cache_table = TRUE, 
    output = "wide") %>% # or tidy
    st_transform(2926)  %>%
    mutate(area_ft = as.numeric(st_area(.)))
## Getting data from the 2013-2017 5-year ACS
colnames(acs5_2018_bg) <- tolower(colnames(acs5_2018_bg))
acs5_2018_bg %>%
    ggplot() +
    geom_sf(aes(fill = medfamincomee), size = .25) +
    scale_fill_viridis_c() + 
    theme_void()

busstop <- st_read(
    file.path(mydatadir, "busstop/busstop.shp"), quiet = TRUE) 
st_crs(busstop) <- 2926
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that
colnames(busstop) <- tolower(colnames(busstop))

print(colnames(busstop))
##  [1] "objectid"   "busstop_id" "zonekey"    "begin_date" "end_date"  
##  [6] "tlink_id"   "ramkey"     "access_lvl" "landing_pd" "trnsfr_pt" 
## [11] "authority"  "awning"     "bearing_cd" "curb"       "curb_ht"   
## [16] "bay"        "bike_rack"  "created_by" "cross_stre" "curb_paint"
## [21] "dt_created" "dt_mapped"  "dt_mod"     "displacemt" "rampositio"
## [26] "ditch"      "ext_surfc"  "ext_width"  "frm_cross"  "frm_intrst"
## [31] "juris"      "reg_fare"   "rfp_dist"   "zip_code"   "i_sgn_anch"
## [36] "i_sgn"      "int_loc"    "link_len"   "pcnt_from"  "t_nd_from" 
## [41] "mod_by"     "news_box"   "non_mt_sgn" "bollards"   "num_shelt" 
## [46] "on_street"  "othr_cov_a" "owner"      "paint_len"  "pk_stp_sur"
## [51] "pullout"    "ret_wall"   "rfa_flag"   "rt_sgn_tp"  "sched_hold"
## [56] "shdr_surf"  "shdr_width" "side"       "side_cross" "side_on"   
## [61] "swlk_width" "sgn_mt_dir" "sgn_pst_an" "sgn_pst_tp" "spc_sgn_tp"
## [66] "length"     "status"     "stop_type"  "address"    "add_commnt"
## [71] "strp_width" "t_signal"   "wlk_surf"   "xcoord"     "ycoord"    
## [76] "xcoord_off" "ycoord_off" "geometry"
busstop <- busstop %>% st_join(acs5_2018_bg)

print(colnames(busstop))
##  [1] "objectid"      "busstop_id"    "zonekey"       "begin_date"   
##  [5] "end_date"      "tlink_id"      "ramkey"        "access_lvl"   
##  [9] "landing_pd"    "trnsfr_pt"     "authority"     "awning"       
## [13] "bearing_cd"    "curb"          "curb_ht"       "bay"          
## [17] "bike_rack"     "created_by"    "cross_stre"    "curb_paint"   
## [21] "dt_created"    "dt_mapped"     "dt_mod"        "displacemt"   
## [25] "rampositio"    "ditch"         "ext_surfc"     "ext_width"    
## [29] "frm_cross"     "frm_intrst"    "juris"         "reg_fare"     
## [33] "rfp_dist"      "zip_code"      "i_sgn_anch"    "i_sgn"        
## [37] "int_loc"       "link_len"      "pcnt_from"     "t_nd_from"    
## [41] "mod_by"        "news_box"      "non_mt_sgn"    "bollards"     
## [45] "num_shelt"     "on_street"     "othr_cov_a"    "owner"        
## [49] "paint_len"     "pk_stp_sur"    "pullout"       "ret_wall"     
## [53] "rfa_flag"      "rt_sgn_tp"     "sched_hold"    "shdr_surf"    
## [57] "shdr_width"    "side"          "side_cross"    "side_on"      
## [61] "swlk_width"    "sgn_mt_dir"    "sgn_pst_an"    "sgn_pst_tp"   
## [65] "spc_sgn_tp"    "length"        "status"        "stop_type"    
## [69] "address"       "add_commnt"    "strp_width"    "t_signal"     
## [73] "wlk_surf"      "xcoord"        "ycoord"        "xcoord_off"   
## [77] "ycoord_off"    "geoid"         "name"          "medfamincomee"
## [81] "medfamincomem" "area_ft"       "geometry"
# tabulate the count of transit stops
nbusstop <- busstop %>% 
    group_by(geoid) %>% 
    summarise(n_busstop = n(), 
              density_ha = n() / min(area_ft) * 107639 , 
              medfamincomee = min(medfamincomee))

nbusstop %>% ggplot(aes(x = medfamincomee, y = density_ha)) + 
    geom_point() + 
    geom_smooth(method = "lm") +
    xlab("block group median family income, ACS-5, 2018") + ylab("transit stop density per ha")
## Warning: Removed 66 rows containing non-finite values (stat_smooth).
## Warning: Removed 66 rows containing missing values (geom_point).

pander(
    summary(
        lm(data = nbusstop, medfamincomee ~ density_ha)))
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 112170 1854 60.49 0
density_ha 89.16 10775 0.008275 0.9934
Fitting linear model: medfamincomee ~ density_ha
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
1198 47100 5.725e-08 -0.0008361

Polygon-on-polygon

nhood <- st_read(
    file.path(mydatadir, "Community_Reporting_Areas.shp"), 
    quiet = TRUE) %>%
    st_transform(2926)
names(nhood) = tolower(names(nhood))

# get the data and project it to match the bus stops, also calculate the area
acs5_2018_trt <- get_acs(
    year = 2018,
    geography = "tract",
    variables = c(n="B06012_001", n_pov="B06012_002"),
    state = "WA",
    county = "King",
    geometry = TRUE,
    moe = 95,
    cache_table = TRUE, 
    output = "wide") %>%
    st_transform(2926) %>%
    mutate(area_ft_tract = as.numeric(st_area(.)))
## Getting data from the 2014-2018 5-year ACS
colnames(acs5_2018_trt) <- tolower(colnames(acs5_2018_trt))
nhood_trt <- st_intersection(x = nhood, acs5_2018_trt) %>% 
    mutate(area_ft_intersect = as.numeric(st_area(.)),
           n_est = ne * as.numeric(st_area(.)) / area_ft_tract, 
           n_est_pov = n_pove * as.numeric(st_area(.)) / area_ft_tract)
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries
nhood_pov <- nhood_trt %>% 
    group_by(gen_alias) %>% 
    summarize(
        neighdist = first(neighdist),
        n = sum(n_est), 
        n_pov = sum(n_est_pov), 
        pct_pov = round(sum(n_est_pov) / sum(n_est) * 100, 1))


nhood_pov %>%
    ggplot() +
    geom_sf(aes(fill = pct_pov), size = .25) +
    scale_fill_viridis_c() + 
    theme_void()

nhood_pov %>%
    ggplot(aes(x = reorder(gen_alias, pct_pov), y=pct_pov)) +
    geom_bar(stat = "identity") +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
    xlab("neighborhood") + ylab("% living under\nfederal poverty level")

st_write(obj = nhood_pov, 
         dsn = mygpkg, 
         layer = "nhood_pov", 
         quiet = TRUE, 
         update = TRUE, 
         delete_layer = TRUE)

st_write(obj = 
             st_union(kcfwy_buf_500ft), 
         dsn = mygpkg, 
         layer = "kcfwy_buf_500ft_union", 
         quiet = TRUE, update = TRUE)
# summarize == union
districts <- nhood_pov %>%
    group_by(neighdist) %>%
    summarise(
        n = sum(n), 
        n_pov = sum(n_pov),
        pct_pov = round(sum(n_pov) / sum(n) * 100, 1))

# save
st_write(obj = districts, dsn = mygpkg, layer = "districts", quiet = TRUE, update = TRUE, delete_layer = TRUE)

Leaflet map

mapview(nhood_pov, zcol = "pct_pov", legend = TRUE)
# CRS
nhood_pov_4326 <- nhood_pov %>% st_transform(4326)

# make a label with the count of persons, persons below poverty, and % poverty
nhood_pov_4326$mylab <- sprintf("n_pov=%s<br>n=%s<br>%s<br>", 
                                round(nhood_pov_4326$n_pov, 0), 
                                round(nhood_pov_4326$n, 0), 
                                paste("pov=", nhood_pov_4326$pct_pov, "%", sep=""))

l <- leaflet(data = nhood_pov_4326) %>%
    addPolygons(popup = ~mylab, weight = 2) %>%
    addTiles()
l
# CRS
nhood_pov_4326 <- nhood_pov %>% st_transform(4326)

mypalette <- colorQuantile(palette = "viridis", domain = nhood_pov_4326$pct_pov, n = 4)

# make a label with the count of persons, persons below poverty, and % poverty
nhood_pov_4326$mylab <- sprintf("%s<br>%s<br>", 
                                nhood_pov_4326$gen_alias, 
                                paste("pov=", nhood_pov_4326$pct_pov, "%", sep=""))

l <- leaflet(data = nhood_pov_4326) %>%
    addPolygons(popup = ~mylab, 
                weight = 2, 
                fillColor = ~mypalette(pct_pov),
                opacity = 0.7) %>%
    addTiles() %>%
    addLegend(pal=mypalette, values=~pct_pov, opacity=0.7, title = "% below poverty", position = "bottomleft" )

l
# CRS
nhood_pov_4326 <- nhood_pov %>% st_transform(4326)

# hospitals
hospitals <- st_read(file.path(mydatadir, "medical_facilities", "medical_facilities.shp"), quiet = TRUE) %>% st_transform(4326)

mypalette <- colorQuantile(palette = "viridis", domain = nhood_pov_4326$pct_pov, n = 4)

# make a label with the count of persons, persons below poverty, and % poverty
nhood_pov_4326$mylab <- sprintf("%s<br>%s<br>", 
                                nhood_pov_4326$gen_alias, 
                                paste("pov=", nhood_pov_4326$pct_pov, "%", sep=""))

l <- leaflet(data = nhood_pov_4326) %>%
    addPolygons(popup = ~mylab, 
                weight = 2, 
                fillColor = ~mypalette(pct_pov),
                opacity = 0.7) %>%
    addTiles() %>%
    addLegend(pal=mypalette, values=~pct_pov, opacity=0.7, title = "% below poverty", position = "bottomleft" )

l <- addCircleMarkers(map = l, 
                     data = hospitals,
                     radius = 5, 
                     weight = 1,
                     opacity = 0.9,
                     fillOpacity = 0.5,
                     label = ~ABB_NAME)

l